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Abstract 

A calculation of neutrino decoupling in the early Universe, including full 
Fermi-Dirac statistics and electron mass dependence in the weak reaction 
rates, is presented. We find that after decoupling, the electron neutrinos 
contribute 0.83% more to the relativistic energy density than in the stan- 
dard scenario, where neutrinos are assumed not to share the heating from 
annihilation. The corresponding number for muon and tau neutrinos is 
0.41%. This has the consequence of modifying the primordial 4 He abundance 
by Ay = +1.0 X 10~ 4 , and the cosmological mass limit on light neutrinos by 
0.2-0.5 eV. 
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I. INTRODUCTION 



The phenomenon of non-equilibrium thermodynamics is of fundamental interest to the 
understanding of the early evolution of the Universe. One example of this non-equilibrium 
behaviour is the decoupling of neutrinos from the electromagnetic plasma, at a temperature 
of a few MeV. Usually, it is assumed that the neutrinos decouple before the temperature 
gets below the mass of electrons and positrons. Thus they do not share the entropy transfer 
from e^ 1 to photons. This means that, although the neutrino distribution is still described by 
an equilibrium distribution after decoupling [1], the temperature is only (4/11) 1 / 3 % 0.714 
times the photon temperature. 

However, the fact that the decoupling temperature is of the same order as the electron 
mass, means that neutrinos must, to a small degree, share the entropy transfer, leading to 
what is known as neutrino heating. This observation has led to several re-investigations of 
neutrino decoupling [2-6]. All of these find a change in the final energy density of neutrinos 
relative to the standard case of total decoupling of around 1%. One problem is, however, that 
they all assume that Boltzmann statistics is adequate for treating this problem. Furthermore 
they assume, that the electron mass is negligible in the weak interaction matrix elements. An 
approximation that is definitely valid whenever the temperature is very large, but possibly 
not when it is comparable to m e . It is therefore of relevance to investigate whether the 
results obtained by using these approximations are correct. 

Since the neutrino distributions are important both in Big Bang nucleosynthesis calcula- 
tions and structure formation models, it is of considerable interest to see if neutrino heating 
has any consequences for either. 

II. FUNDAMENTAL EQUATIONS 

The fundamental equation governing the evolution of all particle abundances is the Boltz- 
mann equation 
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Cl[/] = ECco1i[/], 



(1) 



where Cl is the Liouville operator and C co ii is the sum of all possible collisional interac- 
tions. In the following we will assume that all distributions are homogeneous and isotropic, 
that is /(x,p,t) = f(p,t). In that case, the Liouville operator is [7] 



In the case of standard massless Dirac neutrinos, there are several contributions to C co ii 
[2,8]. To first order in the weak coupling constant all interactions are 2 particle scatterings 
and annihilations. All of these reactions and the corresponding matrix elements have been 
summarized in Tables I and II. Higher order contributions as well as electromagnetic cor- 
rections [9] to these lowest order terms are neglected in the following. Since there are only 
2-particle interactions like 1 + 2 — > 3 + 4, C co ii can be written as 



where A(/i, f 2 , f 3 , U) = (1 - /i)(l - f2)hh ~ (1 ~ /s)(l - fijhfi is the phase space factor, 
including Pauli blocking of the final states, and d 3 p = d 3 p/((2ir) 3 2E). S is a symmetrization 
factor of 1/2! for each pair of identical particles in initial or final states [10], and \M\ 2 is 
the weak interaction matrix element squared, summed over the spin states of all particles 
except the one under scrutiny, pi is the four-momentum of particle i. Note that the before 
mentioned assumptions of Boltzmann statistics and zero electron mass simplify the collision 
integral tremendously, but it is in fact possible to evaluate most of it analytically without 
them. 

Assuming standard weak interactions as described in Tables I and II, it is possible to 
analytically integrate this 9 dimensional integral down to two dimensions (Appendix). This 
makes the integral much easier to evaluate numerically. 

The l.h.s. of the Boltzmann equation contains the factor R~ 1 dR/dt, which cannot be 
immediately evaluated. To do this we have to consider entropy conservation. Let g* be the 






(3) 



S\M\ 2 12 ^ 34 6 4 (p 1 +p 2 -p 3 - P4 )(2Tr) 4 , 
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effective number of relativistic degrees of freedom [1]. If g* is constant in time we obtain 
the standard result of R(t) oc is, but if this is not the case (e.g. during e^ 1 annihilation), 
then this is no longer true. Instead we write R(t) oc f(t), where f(t) can be found using 
the fact that entropy is conserved in the early Universe. In the case of totally decoupled 
neutrinos, f(t) can be calculated without too much trouble. We have used this value for f(t) 
in our calculations as the error induced is small [11]. The photon temperature is calculated 
at each time step by assuming f~ 1 df/dt <C \t~^ ', which is, at all times, an approximation 
better than (9(10~ 3 ). Then the Friedmann equation is 

\dt Rj ~ 4 3 ' [ ) 

and knowing p, p v we can calculate p 1 by assuming that e^ 1 are kept in complete thermody- 
namic equilibrium by electromagnetic reactions with the photon gas. We are thus able to 
follow all relevant thermodynamic quantities as functions of time. 



III. NUMERICAL RESULTS 

To compare with previous results [2,3], we have solved the Boltzmann equation numeri- 
cally for several different cases: 1. m e ^ and FD statistics; 2. m e = and FD statistics; 
3. m e ^ and MB statistics; 4. m e = and MB statistics. In all cases we assume zero 
chemical potentials for the particles involved. 

Case 4 is the approximation used in previous studies. Case 1 is the correct approach. 
The other cases test the errors introduced by the different approximations. For each of 
these scenarios, the energy density deviation, defined as bp v jp v = (p v — p Vo )l p v , has been 
calculated as a function of the photon temperature, T 7 . p Vo is the energy density in a 
neutrino species that has decoupled long before e^ 1 annihilation. For electron neutrinos, the 
result is shown in Fig. 1, whereas Fig. 2 shows the result for muon and tau neutrinos. In 
the limit of MB statistics and zero electron mass (case 4), we reproduce the results of Ref. 
[2] to within 10-15%. This is quite reassuring, as we use a completely different procedure 
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for solving the Boltzmann equations (Dodelson and Turner [2] use a first order perturbation 
expansion). It is seen that, using the correct FD statistics and electron mass (case 1), the 
end result of bp v j p v is around 0.83% for v e and 0.41% for v^. Generally we see that the 
deviation is smaller if FD statistics is used. This is because the phase space factors, A, in 
the collision integral get smaller. Using the correct treatment of electron mass also lessens 
the deviation, due to the changes in matrix element structure. 
Fig. 3 shows the effective neutrino temperature, defined as [2] 

T ' ff(p) s iog(i//(„)-ir (6) 

for FD statistics with and without electron mass. Fig. 4 shows the same, but for MB 
statistics (T e fj = — p/logf). The reason for the offset between the two types of statistics 
is that, for FD statistics, the final neutrino temperature for totally decoupled neutrinos is 
(4/H)i/3 ^ q_ 714 ^ w h ereas f or MB statistics it is (1/3) 1/3 ~ 0.693. We see that, regardless 
of the statistics used, the effective temperature rises with momentum for medium and high 
momentum states. As noted in Refs. [2,3], this is not surprising, because the weak cross 
sections are much larger for large momenta. The shape of the effective temperature curve is 
the same in all four cases, but the actual numerical values are slightly different, being higher 
if the electron mass is neglected in C co u[/]. This effect reflects that the neutrino energy 
density is slightly higher for massless electrons, as seen in Fig. 1. For the lowest momentum 
states the effective temperature rises again. This may seem surprising, but it has to do with 
the integrated matrix elements (the functions called F(pi,p 2 ,P3) in the Appendix). These 
functions contain factors that grow large for very small momenta, leading to a stronger 
interaction of very low momentum states, both with e^ 1 and with high momentum neutrinos. 
This means that the lowest momentum states are actually kept closer to thermal equilibrium 
than the medium momentum states. 
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IV. DISCUSSION AND CONCLUSION 



We have calculated neutrino decoupling in the early Universe by solving the Boltzmann 
equations. Our method differs from previous studies in that it makes no approximations 
for the distribution functions and the electron mass. Although our results do not deviate 
dramatically from these previous calculations, we show that the inclusion of FD statistics 
and non-zero electron mass reduce the effect of neutrino heating on the relativistic neutrino 
energy density by almost 50% for v e compared to the results of Dodelson & Turner [2]. 

We find that the neutrino distribution after decoupling is non-thermal at the 1% level. 
This is also in relatively good agreement with previous results [2,3]. The effective temper- 
ature grows with momentum for medium and large momenta, because of the momentum 
dependence of the weak interactions. The very low momentum states, however, interact 
more strongly, leading to a rise in the effective temperature for very small values of p/T. 

The slight heating of neutrinos relative to the standard scenario has consequences for 
Big Bang nucleosynthesis. We have changed the nucleosynthesis code of Kawano [12] to 
include the effect of neutrino heating. As discussed in Ref. [13], there are several effects that 
combine to change the nucleosynthesis scenario. First of all, the energy density in neutrinos 
is changed. However, this has the consequence of lowering the energy density in photons 
and e ± , because of energy conservation. The result is, if we use FD statistics and non-zero 
electron mass, that the photon temperature goes up because of e^ 1 annihilations by a factor 
of 1.3998 instead of the usual 1.4010 (a change of —0.09%). Furthermore, the weak rates 
are changed, because of the different abundances of e^ 1 and v e . With the relevant changes 
to the code, we obtain a change in the primordial 4 He abundance, Y, of AY = +1.0 X 10~ 4 . 
This is in good agreement with Ref. [13]. They calculate a change of +1.5 X 10~ 4 , but with a 
50% larger neutrino heating. Other authors have found similar values [3,4]. A change in 4 He 
abundance of = +1.0 X 10~ 4 is much below current observational accuracies. The systematic 
uncertainty in the primordial 4 He abundance is estimated to be as large as AY^ ys = ±0.015 
[14], or more than a factor of 100 larger than the change induced by neutrino heating. 
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Another consequence of neutrino heating is that it increases the number density of neu- 
trinos. If one of the neutrino species has a mass, it will therefore contribute slightly more 
than usually expected to f2, the cosmological density parameter. In general, it is possible to 
put a mass limit on a light neutrino that is completely decoupled long before e^ 1 annihilation 
[15]. The mass density of such a species today is 

p v = n v m v . (6) 

Expressing this in terms of the photon density, ra 7 = 2£(3)T^ /ir 2 , one gets 

P, = m^C(3)^. (7) 

If the neutrinos decouple long before e^ 1 annihilation, this can be translated into the normal 
textbook relation 

2 = g, m v 

2 93.03eV ^ ; 

using a present photon temperature of 2.736 K. h is the dimensionless Hubble constant, and 

g v = 2 for one flavor of neutrino and antineutrino. Since observations demand that Q, v < 1, 

we have a mass limit on any given light neutrino. Because of neutrino heating, this limit is 

changed by a small amount. Using FD statistics and non-zero m e , the final number density 

of neutrinos after decoupling deviates from the standard case by bn v jn v ~ 0.52% for v e and 

bn v jn v ~ 0.25% for v^. For electron neutrinos this changes Eq. (8) to 

2 92.55eV ^ ; 

whereas for muon or tau neutrinos we find 

n„ h 2 = 9 ^ T m ^ T . (10) 

2 92.80eV ^ ; 

As this is a very small change to the standard value, neutrino heating does not alter the 
usual conclusion, that electron neutrinos cannot contribute more than about 0.25 (using the 
current experimental upper limit to its mass, m Ve < 7 eV [16]) to $7, whereas muon and tau 
neutrinos can. 



Thus, it is safe to ignore neutrino heating when doing nucleosynthesis and structure 
formation calculations. 
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This appendix shows how to reduce the integral in Eq. (3) from 9 to 2 dimensions. Except 
for the inclusion of mass terms, most of what can be found in this appendix was originally 
developed by Yueh & Buchler [17] for use in supernova calculations. 

The collision integrals all have the form of Eq. (3). The innermost integral can be 
rewritten using the equality 



APPENDIX: COLLISION INTEGRALS 




(Al) 



The integral over d 4 



P4 is now done using the ^-function. Hereafter p 4 is fixed as 



vl = Pi + vl + vl + 2(pi • p 2 ~ Pi ■ P3 ~ P2 ■ Pa). 



(A2) 



Now, the following angles are introduced 



Pi P2 



(A3) 



cos a 



cos 8 



P\P2 
Pi P3 



(A4) 



PlP3 
P2 P3 



(A5) 



cos a 



P2P3 



cos a cos 8 + sin a sin 8 cos (3. 



Thus we obtain 
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d 3 p 2 = p\dp 2 d cos ad/3 (A6) 
d 3 p 3 = p^dp 3 dcos Qd[i. (A7) 



The integration over d/3 is carried out using the ^-function. We use that 



Ji 2 



An important relation for ^-functions is 







i 








d/3 





(A8) 



(A9) 



where /3{ are the roots of f(/3) = 0. The derivative is evaluated to 



d/3 

sin/3i is found as ±(1 — cos 2 /^) 1 / 2 , where 



2p2?3 sin a sin 8 sin/3. (A10) 



2E 2 E 3 -2p 2 p 3 cosacos6 - Q 

cos/3, = : — (All, 

Zp 2 p 3 sin a sin 

(2E 1 E 2 — 2p x p 2 cos a) 



2p 2 ?3 sin a sin # 
(2E 1 E 3 -2 PlP3 cos 6) 



2p 2 p 3 sin a sin 9 

and we have introduced Q = ra\ + ra\ + m?, — m 2 . The equation for sin/3i will have two 
solutions, one in the interval [0,7r], and one in the interval [7r,27r]. Since the r.h.s. of Eq. 
(A9) is symmetric in j3 we can multiply by two and integrate over [0,7r]. The limits on the 
integral dcos a come from demanding that cos 2 (/3) < 1. But this means that 

[2p 2 p 3 sin a sin 6 sin /3) 2 > 0. (A12) 



Notice, that this is the same as demanding that 



Finally we end up with 



2 



> 0. (A13) 
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d/3 


P=Pi 



df3 



(A14) 



P=Pij 



The derivative can be written as 



d(3 



\/ a cos 2 a + b cos a + c, 



(A15) 



/3=/3i 



where 



a = p 2 {~4:K + 8e) 
6 = p 2 (pi -e/ Pl )(87 + 4Q + 8e) 
c = -4 7 2 - 4 7 g - g 2 - 8 7 e - 4Qe - 4e 2 
+4p 2 p 2 (l-cos 2 ^), 



(A16) 
(A17) 
(A18) 



and we have introduced the following parameters in order to limit the amount of space 
required to write out the formulae 



7 = E\E 2 — EiE 3 — E 2 E 3 
e = P1P3 cos 8 



2 | 2 
K = Pi +P 3 - 



(A19) 
(A20) 
(A21) 



Notice that Eqs. (A16), (A17) and (A18) reduce to those of Yueh & Buchler [17] in the 
limit of zero masses. 

Now, any one of the possible matrix elements only include products of 4-momenta. All 
the possible products of these momenta are calculated below: 



Pi ■ p 2 = E 1 E 2 — p x p 2 cos a 
Pi ■ p 3 = E 1 E 3 - p x pz cos 6 
p 1 ■ p 4 = ml + {E X E 2 - p x p 2 cos 

-(E 1 E 3 - p x pz cos 6) 
p 2 - p 3 = [E X E 2 - p x p 2 cos a) 



a 



(A22) 
(A23) 
(A24) 

(A25) 
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CEi£ 3 -PiP3 cos 0) + Q/2 



p 2 ■ P4 = (E 1 E 3 - pip 3 cos 8) + ml 



Q/2 



(A26) 



p 3 ■ p 4 = (E 1 E 2 - pip 2 cos a) - ml + Q/2. 



(A27) 



Because all the above products are analytically integrable over dcos a, the integrals over 
this parameter can now be carried out by use of the fundamental relation 



The step function comes from demanding that there be a real integration interval. This 



interval of [—1, 1], because the step function singles out the physical situations, where dcos a 
cannot be outside this interval. Integrating over d\i is no problem, as there is no dependency 
on this parameter. The final integration that can be done is the one over dcos 8. Any one 
of the possible products of momenta is analytically integrable over dcos 8. The solution of 
b 2 — 4:ac = gives the integration interval. The solutions are 



If there is to be a real integration interval, both of these solutions must be real. The funda- 
mental integration interval is of course [—1, 1], but the real limits are a = sup[— 1, cos #mi n ] 
and /3 = inf [1, cos # max ]. There can only be a real integration interval if both a and /3 are 
real numbers and j3 > a. 

Note that there are seemingly two places where divergences occur in the collision integral. 
The first one is for p\ = p 3 , cos 8 = 1, but in this case there is no problem because although 
the integrand becomes infinite, the integral is finite [17]. The second place is for p\ — > 0. 
To see what happens here, we have to go back to the fundamental integral, Eq. (A28). If 
Pi — > then a = —^p\p\, but then there can be no divergence, as p\ no longer appears in 
any denominator. Clearly there is still the possibility of divergence if m\ — > also, because 




(A28) 



actually also insures that the roots of ax 2 -\-bx-\-c are not outside the fundamental integration 




(A29) 



±2p 2 (2 7 + p\ + p 2 2 + pi + Q)0/(2pip 3 ). 
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of the term. Fortunately it turns out that b 2 — ^ac = if mi,pi = 0, so that the 

rate becomes equal to in this CclSGj clS it ought to be. Finally the collision integral can be 
written as 

2 1 [p 2 2 d P2P ld P3 

Ccoll[/] -J^y W J ~^^ S X ( A3 °) 

A(fi,f2,f3,h)F( Pl , P2 , P3 )e(A), 

where A is the parameter space allowed (that is, the space defined by requiring a,/3 real and 
j3 > a). F is the matrix element integrated over dcos a and dcos 6, 

F{ Pl , P2 , P z) = J \M\ 2 d cos ad cos 8 (A31) 
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FIGURES 

FIG. 1. The evolution oi 6p v /p v for electron neutrinos. The solid curve corresponds to case 1, 
the dashed to case 2, the dotted to case 3, and the dot-dashed to case 4. 



FIG. 2. The evolution of 6p v / p v for muon and tau neutrinos. The curve labels are as in Fig. 1. 

FIG. 3. The effective neutrino temperature after complete annihilation, using FD statistics. 
The solid line is for v e and case 1, the dashed is for v^ T also case 1. The dotted is for v e and case 
2, the dot-dashed for v^ T and case 2. 

FIG. 4. The effective neutrino temperature after complete annihilation, using MB statistics. 
The solid line is for v e and case 3, the dashed is for v^ T also case 3. The dotted is for v e and case 
4, the dot-dashed for v^ T and case 4. 
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TABLES 

TABLE I. Possible electron neutrino processes and the corresponding matrix elements. 1 is 
defined equal to the particle for which we calculate C co ii. 2 is the other incoming particle. 3 
is defined as either particle 1 going out (scattering) or the outgoing lepton (annihilation). 4 is 
defined as either particle 2 going out (scattering) or the outgoing antilepton (annihilation). We 
have introduced the quantities Cy = \ + 2 sin 2 6w and Ca = \, where sin 2 6w ~ 0.23. pi is the 
four-momentum of particle i. 
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TABLE II. Possible fi,T neutrino processes and matrix elements. The definition of particle 

numbers is the same as in TABLE I. 

Process S \M\ 2 

v l + v l ^v l + v l 128G|<(pi • Pi)(p 2 ■ P3) 

Vi + Vj -> Vi + Vj 32G 2 F (p 1 ■ Pi)(p 2 ■ P3) 

f * + f * -> f * + v t 32G 2 F (p 1 ■ p 2 )(p3 ■ Pi) 

Vi + Vj -> Vi + Vj 32G 2 F (p 1 ■ p 2 )(p3 ■ Pi) 

Vi + v e ^> Vi + v e 32G 2 F (p 1 ■ Pi)(p 2 ■ P3) 

Vi + v e -> Vi + v e 32G 2 F (p 1 ■ p 2 )(p3 ■ Pi) 

Vi + e~ ^Vi + e~ 32G 2 F [(C V + C A - 2) 2 ( Pl ■ p 2 )(p 3 ■ p 4 )+ 

{C A ~ C v ) 2 {p\ ■ Pi)(p 2 ■ P3)- 
[(C v -l) 2 -(C A -l) 2 ]m 2 e ( Pl .p 3 )] 
Vi + e+ ^Vi + e+ 32G 2 F [{C V + C A - 2) 2 ( Pl ■ p 4 )(p 2 ■ p 3 )+ 

{C A - C v ) 2 {p\ ■ p 2 )(P3 ■ Pi)~ 
[(C v -l) 2 -(C A -l) 2 ]m 2 e ( Pl .p 3 )] 
Vi + Vi -> e~ + e+ 32G 2 F [{C V + C A - 2) 2 ( Pl ■ p 4 )(p 2 ■ p 3 )+ 

(C A -C v ) 2 (pi-P3)(p 2 -Pi)+ 
[(C v -l) 2 -(C A -l) 2 ]m 2 e ( Pl .p 2 )] 
Vi + Vi -> v e + v e 32G 2 F (p 1 ■ Pi)(p 2 ■ p 3 ) 

32G 2 F { Pl ■ p 4 ){p 2 ■ p 3 ) 
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